function diff = delta_m(m,param)

    % delta(m) in the replacement region

    if param.zeta==0
        diff=param.sigma^2/2/(1-param.alpha).*1./(log(m)+param.B);
    else   
        diff=(param.zeta+param.s.*param.lambda).*q_n(param.m_h,param.psi,param)./q_r(m,param.psi,param)-param.zeta;
    end

    if param.s == 0
        diff=zeros(size(m));
    end

end